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Abstract 



< 

^S^ In this paper, we study the stabihty of various difference approximations of 

C^ the Euler-Korteweg equations. This system of evolution PDEs is a classical 

H isentropic Euler system perturbed by a dispersive (third order) term. The 

'— ' Euler equations are discretized with a classical scheme (e.g. Roe, Rusanov or 

,__! Lax-Friedrichs scheme) whereas the dispersive term is discretized with centered 

E> finite differences. We first prove that a certain amount of numerical viscosity is 

Ji«^ needed for a difference scheme to be stable in the Von Neumann sense. Then 

QQ we consider the entropy stability of difference approximations. For that pur- 

m pose, we introduce an additional unknown, the gradient of a function of the 

■^ density. The Euler-Korteweg system is transformed into a hyperbolic system 

^^ perturbed by a second order skew symmetric term. We prove entropy stability 

of Lax-Friedrichs type schemes under a suitable Courant-Friedrichs-Levy con- 
dition. We validate our approach numerically on a simple case and then carry 
out numerical simulations of a shallow water system with surface tension which 
models thin films down an incline. In addition, we propose a spatial discretiza- 
C^ tion of the Euler-Korteweg system seen as a Hamiltonian system of evolution 

PDEs. This spatial discretization preserves the Hamiltonian structure and thus 
is naturally entropy conservative. This scheme makes possible the numerical 
simulation of the dispersive shock waves of the Euler Korteweg system. 
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1 Introduction 

This paper is motivated by the numerical simulation of the so-called Euler-Korteweg 
system, which arises in the modelling of capillary fluids: these comprise liquid-vapor 
mixtures (for instance highly pressurized and hot water in nuclear reactors cooling 
system, in which the presence of vapor is actually dramatic), superfluids (Helium near 
absolute zero), or even regular fluids at sufficiently small scales (think of ripples on 
shallow water or other thin films). In one space dimension, the most general form of 
the Euler-Korteweg system we consider is 

dtp + d^ipu) = 0, 
^^■^^ dtipu) + d, {pu^ + P{p)) = d, (pt^{p)d,.p + (p^'ip) - «:(p))^^ 



where p denotes the fiuid density, u the fiuid velocity, P{p) the fiuid pressure and 
K,{p) the capillary coefficient. In quantum hydrodynamics, the capillary coefficient 
is chosen so that pK{p) = constant whereas for classical applications, like thin film 
flows, it is often chosen to be constant. In particular, we will focus, for numerical 
computations, to the case of thin film fiows down an incline modeled by the shallow 
water equations with surface tension and source term 

dth + dx{hu) = 0, 
^^•^^ dt{hu) + d^ {hu^ + P(/i, Ai)) = Ai (ghsin{e) - 3^ + -hd^^,h\ , 

where h is the fiuid height, u the streamwise velocity averaged along the cross stream 
direction. The constant p, z/, a are respectively the density, viscosity and capillarity of 
the fiuid used in the Liu-Gollub experiment |LGj whereas g is the constant of gravity 
and 9 is the slope of the channel. The constant Ai > is arbitrary and P{h, Ai) is a 
smooth pressure term. 



The Euler-Korteweg system (1.1) falls in the class of abstract Hamiltonian systems 
of evolutions PDEs when it is written with variables p, u: 

(1.3) dtU = J iY^Hp]) , 

with U = {p^uY , J = dx^, 

and E denotes the Euler operator 

En[u] = I y + ^'^p^ + ^'(p^^-^ - 9Mp)d.p) ) = ( y + ^p^(p^ ^^p) 

pu j \ pu 
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The pressure P is related to F through the relation pF'{p) —F{p) = P{p). Due to the 
invariance of the equations with respect to spatial translations and time translations, 



the system (1.3) admits, via Noether's theorem, two additional conservation laws 
which are nothing but the conservation of momentum (and the second equation of 
(LTJ) and the conservation of energy: 



(1.4) 



dtipu) + d, {pu^ + P{p)) = d, (pK{p)d,,p + {pk\p) - k(p))^^ 

dt { ^pu^ + ^(p, d^p) \ +da,{ -pu^ + puEp£{p, d^p) + K{p)d^{pu)d^p j = 0. 



As a consequence, when the system (1.1) is set on the real line or with periodic 



boundary conditions, the energy (here the Hamiltonian "H) is conserved. Therefore, 



it is desirable from a numerical point of view that a difference approximation of (1.1) 



or (1.3) preserves the energy or, at least, dissipates energy. In the first case, the 
difference approximation is an "entropy conservative" scheme and in the later case, 
it is an "entropy stable" scheme. 

At this stage, there are two possible strategies to tackle this problem. The first one 



consists in considering the system (1.1) as a dispersive perturbation of the classical 
isentropic Euler equations. This is the point of view adopted e.g. in |LMRj . Here the 
authors construct fully discrete entropy conservative scheme for systems of conser- 
vation laws (hyperbolic or hyperbolic-elliptic) endowed with an entropy-entropy flux 
pair. These difference approximations are second and third order accurate and can in 
turn be used to construct a numerical method for the computation of weak solutions 
containing non- classical regularization-sensitive shock waves. In particular, the au- 
thors considered dissipative/dispersive regularizations that are linear in the entropy 
variables 

dtu{v') + dj{u{v')) = sB^d^^v' + e'^B^d^^^v', < e < 1 

with Bi constant matrices, B2 being positive definite and S3 symmetric. Therefore 
the dispersive terms do not contribute in the energy equation: 



dt / U{u') < 0. 

Jr 

where U is the entropy associated to the system of conservation laws 

dtu + d^J{u) = 0. 
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This situation contrasts with the one met in the Euler-Korteweg system where the 
dispersive terms have a contribution in the energy balance 



5.^Qp«'+i^(p))+'^(p)^=o. 



As a consequence, an entropy conservative or entropy stable scheme for the isentropic 
Euler equations coupled with a centered approximation of dispersive terms may not 
provide an entropy conservative nor entropy stable scheme for the Euler-Korteweg 
system. This issue was considered in |CLj where Euler Kortweg equations are written 
in lagrangian coordinates of mass: by introducing an extended formulation of the sys- 
tem, the authors derived a family of high order and entropy conservative semi-discrete 
schemes. With these high order approximations schemes in hand, the authors then 
computed kinetic relations for Van der Waals fluids. Though, the lagrangian coor- 
dinates of mass can not be used in dimension d with d > 2 and one has to consider 
an alternative extended formulation in Eulerian coordinates: this latter point of view 
will be expanded here, based on the extended formulation found in |BDDl IBDDd] . 

In section [2| we consider the stability of various difference approximations of Euler- 
Korteweg equations in the Von Neumann sense. We shall prove that even at that 
linear level, the Godunov scheme (explicit and implicit in time) is always unstable. 
In this direction, we checked the stability of Lax-Friedrichs type schemes: we show 
that it is stable in the Von Neumann sense under a suitable Courant-Friedrichs-Levy 
(CFL) condition for explicit forward Euler (resp. Runge Kutta) time discretization 
for first order (resp. second order) difference schemes. This analysis provides neces- 
sary conditions of stability for the simulations of the fully nonlinear system. Finally 
we show that the implicit backward Euler time discretization is always stable. 



In section pi we move to the entropy (nonlinear) stability problem. It is a hard prob- 
lem to obtain directly entropy stability from nonlinear difference approximation of 
Euler Korteweg equations since discrete integration by parts and time discretization 
do not commute. Here, we introduce an additional variable w = \f^^{fi)dxP I ^ and 
derive a conservation law for w. In this new formulation, the capillary term appears as 
an anti dissipative term in the system for (m, w) and one can prove the well-posedness 
of the Euler-Korteweg system |BDDj . Moreover, the derivation of the energy estimate 
follows the same line as a classical energy estimate in the isentropic Euler equations. 
In that setting, we show difference approximations made of a Lax-Friedrichs (entropy 
stable) scheme for the hyperbolic part and centered difference for the anti-diffusive 
part is entropy stable under a suitable CFL condition for explicit forward Euler time 
discretization and always stable for implicit backward Euler time discretization. 
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Finally, in section [4], we carry out numerical simulations of shallow water equations 
with surface tension. We first consider thin film flow over a flat bottom and neglect 
source terms so as to compare entropy stability of difference approximations for shal- 
low water in original form and for its new formulation counterpart. The numerical 
simulations clearly show that the discretization of new formulation of shallow water 
equations has better entropy stability property. 

Then, we introduce a difference approximation of the shallow water equations 
written as a Hamiltonian system of evolution PDEs. The semi discretized system is 
Hamiltonian and trivially preserves a discrete Hamiltonian that is consistent with the 
continuous one. Then, one is left with the problem of time discretization: the explicit 
forward Euler is always unstable and one has to consider implicit time discretization 
to obtain an entropy stable scheme: the numerical simulation of this Hamiltonian 
system shows that the dynamical behavior is completely changed in comparison to 
entropy stable schemes. Indeed, this Hamiltonian difference approximation has no 
numerical viscosity, so that one can observe the formation of so called "dispersive 
shock waves" [El lEGK^ lEGS] . Here, the classic hyperbolic shocks are regularized by 
dispersive effects and an oscillatory zone appear and grows with time. We conclude 
this section with numerical simulation of a Liu Gollub experiment |LGj modeled by 
a consistent shallow water model with source term derived in |NVj . 

2 Von Neumann stability of difference schemes 

In this section, we linearize the Euler-Korteweg system about constant states and 
study the Von Neumann-stability of various difference approximations. We have first 
considered two classes of spatial discretisations, namely Godunov and Lax-Friedrichs 
type schemes for the first part of the equations whereas the dispersive term is dis- 
cretized with classical centered difference approximations. We prove that the Go- 
dunov scheme is always unstable both with explicit forward and implicit backward 
Euler time discretization. The Lax-Friedrichs type schemes inherit better stability 
properties: we prove that it is stable under a CFL condition for explicit forward Eu- 
ler time discretization and always stable for implicit backward Euler time discretiza- 
tion. These difference approximations are first order accurate in time and space: we 
have also considered second order accurate schemes, namely MUSCL scheme with a 
Lax-Friedrichs flux. We prove the stability of explicit Runge Kutta (second order 
accurate) time integration under a CFL condition and stability of Crank-Nicolson 
time discretization. 

The linearized Euler-Korteweg equations (in (p, q = pu) variables) are given by 
(2.1) dtp + d^q = 0, dtq + (c^ - u'^)d^p + 2ud^q = ad^^^p, 
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with a = pK,{p). With Riemann invariants, this system is also written as 
(2.2) dtr + {u-c)d^r = -adj:j::r{r + s), dts + {u + ajd^s = adxxx{r + s). 

2.1 Instability of Godunov type schemes 

In this section, we consider Godunov schemes for the hyperbohc part of the equations 
and centered finite difference approximation of the third order partial derivative. In 
order to simplify the analysis, we deal with a diagonalized version of the first order 



part (2.2). Let us introduce the matrices A and B 



A=i ^7' ° ), B= -"" -^ 

u+c r \ ^ ^ 



so that, by setting v = (r, s)^, the system (2.2) reads 

(2.3) dtv + AdxV = BdxxxV- 

We first consider supersonic flows u > c > 0: the Godunov scheme reads 

(2.4) ^ + a"^- ~ "^-^ = i?"^+^ ~ ^"^-^^ + ?"^-^ ~ "^-^ , ^eZ. 

dt 6x 2dx^ 

Then the explicit forward Euler (FE) and implicit backward Euler (BE) read 

(2.5) t;;+i - v] + X.Aiv] - v]_,) = ^B (t;;^^ - 2t;;+i + 2v]_, - v]_,) , 

(2.6) v-^' - vj + X.AivJ^' - vpl) = ^B (t;;+' - 2^+^ + 2vp^ - vp^) , 

where A^ = 5t/5x^ encode the classical CFL conditions for hyperbolic equations 
(fc = 1), parabolic equations {k = 2) and dispersive equations (fc = 3). We prove the 
following proposition. 

Proposition 2.1. Assume u > c > 0, there exist A > and h > such that if 



< A3 < A and < 5x < h then the difference approximations (2.5) and (2.6) of the 



linearized Euler- Kortew eg system are unstable in the Von Neumann sense. 

Proof. Let us start with the explicit forward Euler time discretization: we search 
for solution of (2.5) in the form f" = X"'e~^^^v. Then, {X,v) satisfies the spectral 
problem 



^. _ /I + XM){u - c) - ta-fiOXs -^cT7(e)A 



M2{^,Xi,Xs)v, 
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with S{0 = e*« - 1 and 7(0 = 4 sin(0(l - cos(O). Hence X lies in Sp(M2(C, Ai, A3)) 
and satisfies the following equation: 

(X - 1)2 - 2Ai5(Om(X - 1) - 2zac7(05(0AiA3 + OiXJ) = 0. 

Recall that Ai = X^Sx"^. Next, we expand X with respect to < A3 ^ 1 and 
< Sx <t^ 1: for that purpose, we set X = 1 + X^SxX. Then, one has 

X^ - 2zac7(05(0 = 25(0^ 5x + 0(5a;2), 

where 0{Sx'^) denotes a smooth function of 6x, independent of Aj and of order 0{dx'^). 
Then, we choose ^ = 7r/2: by a classic application of the implicit function theorem, 
one obtains that X G Sp(M2(7r/2, Ai, A3)) expands as 

IT 

X = 1 ± 4v/^ccos(7r/4)e*8A3(5a: + OiX^Sx^). 

As a result, one finds 

|Xp = 1 ± 8^/accos{n/4:) cos(7r/8)A35x + OiXsSx^). 

Then, for A3 and 6x sufficiently small, M2(7r/2, Ai, A3) has two eigenvalues X± such 
that |X_| < 1 < 1X4. 1 and the explicit in time Godunov scheme is unstable in the 
Von Neumann sense. 

The instability of the implicit Godunov scheme then follows easily. Indeed, by search- 
ing for solutions of (2.6) in the form f" = X"'e~^^^v, one obtains the spectral problem 



Xv = M2(e, -Ai, -X^y'v = M3(e, Ai, Xs)v, 

and it is easily seen that for A3 and 6x sufficiently small, M3(7r/2, Ai, A3) has two 
eigenvalues X± so that |X_| < 1 < |X+| and the implicit in time Godunov scheme is 
unstable. D 

Remark 2.2. This instability can be explained by the fact that for the continuous 
system the flow is not always supercritical. Indeed, searching for solutions of the 
continuous system (2.1) in the form {p,q) = e*'^*^^~**-*(p, g), one finds 



s±{k) = u± \i (? + ak'^. 



Then, if\k\ < kc = \/{u'^ — (?)/(^, the flow is supercritical .'0 < s^{k) < s^{k) whereas 
for sufficiently large wavenumbers \k\ > kc, the flow is subcritical: S-{k) < < s^{k). 
As a consequence an upwind strategy necessarily fails: instability is found for wave 
trains with large wave number (or equivalently for small 6x). 
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Though, even in the subcritical case < m < c where the flow remains subcritical 
when the dispersive term is added (s_(A;) < < s+(A;), VA; G M), the Godunov scheme 
is unstable. We prove the following proposition. 

Proposition 2.3. For subcritical flows, < u < c, there exists A > and h > 
such that z/ < A3 < A and < 6x < h then the Godunov space discretization yields 
an unstable scheme in the Von Neumann sense for both explicit forward and implicit 
backward Euler time discretization. 

Proof. The semi discretized Godunov scheme is written as 



; +{u c) ' 
dt Ox 


- b""^^' - 


- 2vj+i + 2vj.i 


- Vj.2 


dsj . x^i - r,-i 
I dl^^^^'^ 6x ) 




2Sx^ 





We use the notations introduced in the proof of proposition 2.1 Let us first consider 



the explicit forward Euler time discretization. In this case, the spectral problem reads 



\ 2CT7(e)A3 l + Ai(5(0(M + c)+z(T7(e)A3 

= M2(e,Ai,A3){). 
Then X lies in Sp(M2(^, Ai, A3)) and satisfies 



(X - 1)2 - XMOic - u)+6{0{u + c)){X - 1) 

+ zaAiA37(0 {mic -u)- 5(0(u + c)) + 0(A?) = 0. 

By setting X = 1 + X^SxX, one obtains 

X^ + lajiO (W)ic -u)- 6{0iu + c)) = Sx{6{^{c - u) + 5(0(n + c))X + 0{6x^). 

This is also written as 

X^ + 2^7(0 («m(1 - cos(O) + csin(O) = 0{Sx). 

Let us choose ^ = 7r/2: then, it is easily proved that X G a{M2{TT/2, Ai, A3)) expands 

as 

n 

X = 1 ± 4v/accos(7r/4)e*8A3(5x + 0{XsSx^). 

Thus, one easily finds that for A3 and dx sufficiently small, there are two eigenvalues 
X± so that |X_| < 1 < |X+| and the explicit Godunov scheme is unstable. By using 
a similar argument, one proves easily that the implicit backward time discretization 
yields an unstable scheme as well. D 
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Remark 2.4. The instability still arises in the subcritical case since the Riemann 
invariants r, s for the first order system does not always represent respectively a sub- 
critical and supercritical mode: it depends again on the wavenumber. 



2.2 Stability of Lax-Friedrichs type schemes 

In this section, we check the stabihty of Lax Friedrichs type scheme in the Von 
Neumann sense: this provides necessary, and in general, practical conditions, for the 
stability of the nonlinear difference approximations. Since there is no upwind strategy. 



we discretized directly (2.1). Let us introduce the matrices A and B 



A 







c2 



u^ 



1 
2m 



B 




a 



so that system (2.1) reads 
(2.7) 



dfV + AdxV = BdxxxV- 



2.2.1 First order accurate schemes 

We first discretize in space by using centered differences for the third order deriva- 
tive and a Lax Friedrichs type flux for the hyperbolic part. The semi-discretized 
approximation then reads 



(2.8) 



dt 



+ ;r^^^ {vj+2 - 2vj+i + 2vj.i - Vj^2) 



_^.^i+i - ^i-i 



+ 



X 



2fe3 



25x 
where A+Vj = Vj+i — Vj and 

gLF{ui,Ur) 



— yVj+i-2vj + Vj^i) + B^' 



Vj+2 - 2f j+1 + 2Vj-i - Vj.2 



2(5x3 



F iui) + F iUr 



Q_ 

2Ai 



{Ur - Ui) 



is a general Lax-Friedrichs flux. The classical Lax Friedrichs scheme is obtained by 
setting Q = 1 , whereas ii Q = \ one finds the modified Lax-Friedrichs scheme [T]. 
Finally, by setting Q = Ai ||y4|| = Ai('U + c), one obtains the Rusanov scheme. 



We search for spatially bounded solutions of (2.8) in the form Vj = e^^^v. This yields 
(2.9) 



dij 
'dt 
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with 



e 



-*5 _ pii 



«(e) = 2z = " '^''(^)' ^(^) = ^ -2 + e^^ = 2(cos(0 - 1), 

7(0 = ^ = 2sin(0 - sin(20 = 2sm(0(l - cos(O). 

Let us now consider various time discretizations. The forward exphcit /backward 



iniphcit Euler time discretizations and ^-scheme for (2.9) are written as 

• (FE) Euler: t)"+i - {)" + {i\ia{^)A - f /3(0ld){)" = ^37(0^^^" 

• (BI) Euler: {)"+i - {)" + {tXia{OA - f /3(0ld){)"+^ = 2A37(05{)"+i 

• ^-Scheme: v^'^^ = {)" - (iXiai^A - ^/3(0ld - ^37(0^ j ((1 - e)t)"+^ + ev"") 

Theses difference approximations are all first order accurate both in space and time. 
In the particular case 6 = 1/2 (Crank-Nicolson time discretization), the difference 
approximation is second order accurate in time. 

We first study the (FE) time disretization. Denote R{Q,m, \u\,c) the function 

R[Q,m,\u\,c) = mm ? 

o<^<2 \^(^2 - x) {\u\ + Vc^ + 2mxy 

Note that 

1 

R (l,m, |m|, c) = ■ 



(\u\ + a/c2 + Am) 
We prove the following proposition. 

Proposition 2.5. The (FE) time discretization yields a stable scheme in the Von 
Neumann sense if and only if < Q < 1 and 

Proof. The (FE) time discretization yields the finite difference scheme: 

{;"+! = (1 - Q(l - cos(0))ld + z sin(OMo(e, Xi, A3)) i)" 

with Mo(^, AijAs) = XiA + 2A3i?(l — cos(^)). The scheme is stable in the Von 
Neumann sense provided that 

Sp (((1 - g(l - cos(0))ld + ^ sin(OMo(e, Ai, A3))) C {A G C; |A| < 1}. 
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By an easy computation, one shows that 



Sp (Mo(e, Ai, A3)) = <; Ai ( n ± ^/c2 + 2a^(l - cos(O) | } = {Ai7±}. 

Let X = 1 — cos(^) : the Lax Friedrichs scheme is stable if 

2 



;i - Qxy + Xix (2 - x) I \u\ + \lc'^ + 2^Y^a; ) < 1. 

A necessary condition for the stabihty of the finite difference scheme is < Q < 1- 
Let us now define 

K[Q,m,\u\,c) = mm -, 

o<x<2 \(2-x) {\u\ + Vc^ + 2mxy 

Then, the finite difference scheme is stable under the assumption < Q < 1 and 

Ai < R{Q,a—, |u|,c). 
Ai 

This completes the proof of the proposition. D 

Corollary 2.6. The difference approximation associated to the classical Lax- Friedrichs 
scheme, Q = 1, is stable in the Von Neumann sense if and only if 



Ai I H + ^/c2 + ^j<l (CFL),. 

The difference approximation associated to the Rusanov scheme, Q = {u-\-c)Xi = pXi 
is stable in the Von Neumann sense if 



Ai^ ^ — = '— < 1- 

\u\ + c 

Remark 2.7. In the case of the classical Friedrichs scheme, we obtained a Courant- 
Friedrichs-Levy (CFL) condition for 6t = 0((5x^). Indeed, one can see heuristically 
that this condition is a necessary condition of stability. Indeed, in Fourier space. 



the wave speeds of (2.7) are written as s{k) = u ± y/c^ + a/c^ . Then, in the large 



wavenumber limit, one has s{k) < C\k\. Heuristically, it is necessary for a numerical 



2 VON NEUMANN STABILITY OF DIFFERENCE SCHEMES 12 



stable to be stable that the domain of dependence of the numerical solution contains 
the domain of dependence of the exact solution. This is also written as s{k) 6t/6x < 1. 
On a spatial grid with stepsize Sx, one has s{k) < C/6x since the largest wavenumber 
is 0{l/6x). As a consequence, one obtains a heuristic CFL condition Cbtjbx^ < 1 
which is precisely the CFL type condition for the Lax-Friedrichs scheme. It is easily 
seen on the Rusanov scheme that this condition is not sufficient to get stability: here 
a necessary condition of stability is 6t = 0{6x^). 



Let us now consider backward implicit Euler time discretization of (|2.9|): 

(1 + Q(l - cos(0))ld - z sin(OMo(e, Ai, A3)) v 
This scheme is stable in the Von Neumann sense if and only if 



\2 , ^2■2,.^ I ,-.i , /^2 , o^-^3, 



;i + g (1 - cos(O))' + xi sin^(e) ( 1^1 + y c^ + 2a^(i - cos(O) 1 > i, 

which, obviously, holds true for all Q > 0. 

We end this section by checking the stability of ^-schemes. We prove the following 
proposition. 

Proposition 2.8. If 9 > ^, the 9-scheme is stable in the Von Neumann sense if and 
only if 

{20-l)Q<l 

with 



{2e-iYxi<R\Q{2e-i),^,\uic 

If Q ^\, the 6 -scheme is always stable. 

Proof. The ^-scheme can be written as -0"+^ = Mi{C,, Ai, A3, 6)v"' and the spectrum of 
A^i(^, Ai,A3,6') is given by 



1 + (1 - 9) 0(1 - cos(O) - i (1 - 9) smK)An± 
As a consequence, the ^-scheme is stable if and only if 



:i + (1 - e)Qil - cos(0))2 + (1 - e)^ sin^(0(Ai7±)2 



< 1. 
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This condition is also written as 

{29 - 1) sin2(0(Ai7±)' < (2 + (1 - 2e)Q{l - cos(0))g(l - cos(O). 
Since |7+| > |7_|, the 6'-scheme is stable if 

(26 - 1) sin2(0(Ai7+)' < (2 + (1 - 2^)Q(1 - cos(0))Q(l - cos(O). 

As a consequence, the ^-scheme is always stable if ^ < 1/2. Let us now consider the 
case 6 > 1/2: then one has the following CFL condition: 

2,2 / 2 - (20 - l)g(l - cos(O(20 - l)g(l - cos(e)) 



{26 - lyxi < 



[1 + cos{o){\u\ + sic' + 2^(1 - ^^<o)y 



The discussion is now similar to the one carried out in the explicit case 6 = 0, namely 



\2 x2 /D / ^/OZ) i\ -^3 



{26 - ly Xi<RiQ {26 - 1) , T^, \u\, c 

This concludes the proof of the proposition. D 

Remark 2.9. Note that for the Crank Nicolson scheme, 6 = 1/2, one can choose 
Q = a classical centered scheme. In this case, this corresponds to the numerical 
schemes used for the practical simulation of thin film flows down an inclined plane in 
the presence of surface tension jiKRSV^ . 

2.2.2 Second order accurate schemes 

Hereafter, we consider second order accurate schemes and consider Lax-Friedrichs 
fluxes . For the time discretization, we consider the Runge Kutta (second order 
accurate) method and the Crank Nicolson method. Assume that the linearized Euler- 
Korteweg equations are written as 

vt + Ad^v = Bdn^xxV- 

We discretize in space by using a MUSCL scheme |Colt IVLj for the first order dif- 
ferential operator without nonlinear monotony correction of the slope (it does not 
operate in the smooth monotone area of the solution), and centered approximation 
of third order differential terms : 



dVj A+gLF{vj-l/2-,Vj-l/2,+) , B o ,o \ 
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where f,-i/2,- = Vj-i + -^, v 



di 



2-5 ^i-i/2,+ - ^i - T ) 



dj is the local increment of v given 



(without monotonicity correction) by the centered formula dj 



^i+i-^j-i 



dvj _ A {vj+2 - Qvj+i + Qvj-i - Vj-2) Q i-Vj+2 + 4'Uj+i - Qvj + Avj^i - ■Uj_2) 
'dt ~ S6x ^2X1 S5x 



(2.10) 



B 



^25^ ^'"^^^ ~ '^^''^^ ^ '^^^'^ ~ ^^'^^ 



with j G Z. We search for solutions of (2.10) in the form Vj{t) = e ^^^v{t), one finds 



with 



and 



dv 



Mo(e, Ai, A3) = y ^ (3 - cos(O) + 2A3i? (1 - cos(O) 



M{C,Q,\iA3 



Q. 



1 - cos(0)'W + t sin(OMo(e, Ai, A3). 



Then, the second order accurate scheme with explicit Runge-Kutta time discretization 

is written as 



-,n+l 



1 



2 1 „'-,« 



Id + M{^, Q, Ai, A3) + -M{^, Q, Ai, A3)' V' 



The eigenvalues of Ai{^, Q, Ai, A3) are given by A-t = — ^(1 — cos(,^))^ — zAi sin(^)7-|- 
with 

3 — cos(^) 



7± 



_^ _ 4aA3(l^cos(0) 
^'^ Ai(3-cos(0) 



In order to simplify notations, we introduce Q = y (1 — cos(^))' and Ai = Ai7± sin(^). 
Then, the second order accurate scheme with Runge Kutta time discretization is 
stable if and only if 



-Q-iAi 
'l-Q-iAi + ^ ^— <1 



or equivalently 



Y<Q-^ + \/2Q-Q' 
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2 + a;\ / _ / 4(7 X^x 



With Q = \Qx and A^ = \\x (2 - x) ( ^^ ) ( u + Ac^ + 



Ai(2 + a:) 
As a result, the exphcit second order accurate scheme (MUSCL) is stable only if 



A?h+J--^ + 7J?^l <4 



' v/4g - g2^2 + g^ _ ^ 



(5x2 (2 + a;)/ - (2-a;)(2 + x)2 

and we have proved the following proposition. 

Proposition 2.10. The second order accurate scheme with MUSCL type discretiza- 
tion in space and Runge Kutta time discretization is stable if and only if < Q < 1 
and 

for all X G [0,2]. 

Corollary 2.11. The classical Lax-Friedrichs scheme with MUSCL space discretiza- 
tion is stable in the Von Neumann sense if 



V4g 


-Q' 


'■x^ + Qx- 


Q-'x-' 
4 




(2- 


x){2 + xY 





AiU+a/c2 + ^1 <i 



When 5x is sufficiently small one finds the following condition, : 



6t [V^±RY7'^VS ( r- 



_ ^ I V l^ + cj j ' "- 
5x3 

for the Rusanov scheme with MUSCL space discretization 

Remark 2.12. For the classical Lax Friedrichs scheme, the CFL stability condition 
has the same nature than in the case of first order accurate schemes whereas we get an 
improved CFL condition for the Rusanov scheme in comparison to first order accurate 
schemes. 

We finish this section by checking the stability of an other second order scheme, 
namely the Crank Nicolson scheme (^-scheme with 9 = 1/2). We prove the following 
proposition. 
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Proposition 2.13. The difference approximation with Crank Nicolson time discretiza- 
tion and second order accurate in space (MUSCL with Lax-Friedrichs fluxes) is stable 
for all Q > 0. 



Proof. This Crank-Nicolson scheme for (2.10) reads 

1 



-.n+l 



1 



:A^(e,Q,Ai,A3 



-1 



1 



M{i,QAiA?. 



It is stable if and only if 



l_2(l_cos(0)' + ^^sin(07± 



1 + ^(1 



cos(0)2-^f sin(07± 



< 1 



It is easily seen that this condition is equivalent to 

(4 - g(l - cos(O)')' + (2Ai7±)' < (4 + g(l - cosiOff + (2Ai7±)', 

which obviously holds true for any Q > 0, and concludes the proof of the proposition. 

D 



3 Entropy stability of difference approximations 

In this section, we study the entropy stability of difference approximations for Euler- 



Korteweg equations (1.1). Recall that (p, m) solution of (1.1) satisfies the energy 
estimate 



u 



dt P^ + F{p) + K{p) 



{d..pY 



0, 



with T = M./LZ for any L > 0. The surface tension plays a significant role in the 
energy estimate and the previous section illustrates that it is a non trivial task to 
obtain a numerical scheme which conserves or, at least, dissipates the energy, even at 
the linearized level. 

In this section, we introduce a new unknown w = \/K.{p)dxp/ y/p and derive an 
evolution equation for w. The system of evolution PDEs for (p, pu, pw) is made of a 
first order hyperbolic part perturbed by a second order anti dissipative term. This 
latter term is discretized by centered finite differences. We show that any entropy 
dissipative schemes for the hyperbolic part (in the sense defined by Tadmor in [T]). 
provides an entropy dissipative scheme for the "augmented" Euler-Korteweg system. 
Then we consider fully discrete schemes and show that (FE) time discretization is 
entropy stable under a suitable CFL condition which is consistent with the one derived 
in the previous section whereas (BE) time discretization is always stable. 
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3.1 "Schrodinger" formulation of the Euler-Korteweg system 

We start from the Euler-Korteweg system written in (p, pu) variables: by setting 
■"^ = \I~^Ap)^xPI \fp-, one finds the augmented Euler Korteweg system 



(3.1) 
(3.2) 
(3.3) 

with \i{p) 



dt{pu) + 9^.(pM^ + P{p)) = d^ {p{p)d^w) , 
dtipw) + d^{puw) = -d^ {p{p)d^u) 



P 



3/2, 



k(p). The equation (3.3[) is derived from the mass conservation 



law (3.1). In order to simplify notations, we introduce v = {p, pu, pw) and f{v) 



{pu, p'ir + P{p), puw)^: then the system (3.1 3.3) reads 



(3.4) dtV + dJiv) = d,iBip)d^z), 

where z = p~^v and B{p) denotes the skew-symmetric matrix 



B{p) = ( Pip) 
-Pip) 



Note that the first order part of (3.4) (-B(p) = 0) admits an entropy-entropy flux pair 
(f/, G) with 



Uiv) = Fip)+p 



u^ + w"^ 



Giv) = uiUiv) + Pip)) 



whereas the augmented Euler-Korteweg system (3.4) admits an additional conserva- 
tion law 

dtUiv) + d^Giv) = d^ ipip)iud^w - wd^^u)) . 



We consider difference approximations of (3.4) in the conservative form 






fe2 



Bipj+i) izj+i - Zj) - 5(p._i) izj - Zj^i)) . 



Following the terminology used in [T], we enquire when the difference schemes ( |3.5 ) 
are entropy stable in the sense that there exists a numerical flux Q^.i, which takes 

-'"'"2 



also into account the right hand side of (3.5), so that 



(3.6) 



dt "' dx 



< 0. 



The difference approximation (3.5) is entropy conservative if the inequality in (3.6) 



is an equality. Note that any entropy-stable scheme satisfies the entropy inequal- 



ity of the original system (1.1) in a weaker sense since Wjit) is an approximation of 
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yf^^ipjdxp/y/p at the point Xj = j 5x. In the last section of the paper, we will use 
the Hamiltonian structure of (1.1) to obtain an entropy conservative scheme (for the 
semi discrete problem). 

In what follows, we prove the following proposition. 
Proposition 3.1. Let us consider the finite difference scheme 

d ,., fj+i-fi- 



(3.7) 



dt 



v,it) + 



'3-\ 



dx 



0, 



difference scheme (3.5) is entropy stable. 



which is a semi discretization of (3.4) with B = and is entropy stable. Then the 



Proof. The difference approximation (3.7) is entropy stable, then there exists a con- 
sistent entropy flux J^j+i so that 



U.{v,f^^^^ ^- 



J", 



i+ 



-F. 1 



J-; 



6x 



6x 



+ n 



•J' 



with TZj > 0. Let us now consider the difference scheme (3.5) and multiply this 
equation by U^ivjY: one finds 



d 
dt 



U{v,) + 



5x 



+ ^. = ^^ (5(P,+|) (^.+1 - ^.) - 5(P,-|) (^. - 



Zj-l 



:=/C 



V 



We focus on the "capillary" term Kj: it is written as 

Sx'^JCj =Uj ifij^i{wj^i — Wj) — fij_i{wj — Wj-i) ) 

-Wj (^fij^{Uj+i -Uj)- f^j,i{Uj - Mj_i)j 

(3.8) 
Then, we introduce the entropy flux Gj^i'- 

yj+h = •^:m - ^'^H ji ■ 

This numerical entropy flux is clearly consistent with the continuous entropy flux of 
the augmented Euler Korteweg system 



^(p, M, w, d^u, d^w) = u {U{v) + P{p)) - p{p) {ud^w - wd^u) 
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provided the entropy flux J^j^i is consistent with the entropy flux G{v) of the hy- 



perbohc part of (3.1 3.3). Moreover, we have the following semi discrete entropy 
estimate 



d 



U{v, 



^j+l 



dt ' -" 5x 

This completes the proof of the proposition. 



-7^J < 0. 



D 



By applying the proposition (3.1) and considering the analysis of various entropy 



stable schemes found in [T], one finds that many of the classical three points (first or- 
der) schemes (Rusanov scheme, Lax Friedrichs scheme, Harten-Lax-van Leer scheme) 
provide natural entropy stable schemes for the augmented Euler-Korteweg system 



((3.1 )-(3.3)). For apphcation purposes, we check the entropy stability oi fully discrete 



schemes. 



3.2 Entropy stability of fully-discrete schemes 

In this section, we restrict our discussion to the backward implicit and the forward 
explicit Euler time discretization which read, respectively: 



(3.9) 



(3.10) 



,..1 



^" + ^^(^-^) 



^AB{pZ\){ 






,..1 



•^ 2 



We first prove the entropy stability of the implicit backward Euler time discretization. 
Proposition 3.2. Let us consider the semi discretized scheme 



(3.11) 



dt'^^^ dx 



0, 



which is an entropy stable approximation of (3.4) with B = 0. Then the associated 



(BE) time discretization (3.9) is always entropy stable: there exists Q\i so that 



(3.12) 



6t 



u{v^+') - u{v]) + -{g;^, - gj_,j < o, vj, vn. 
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Proof. The entropy f/ is a convex function of the variables v = {p,u,w): then one 
has 

(3.13) u{vj+') < u{v]) + u,{v;+Y{v]+' - v:;). 

The semi-discrete scheme ( |3.11[ ) is entropy stable so that 

with TZ^ > 0. Moreover, one has 

(3.14) = p^ll {u]^'w;tl - u]tlw;+') - /x;+l {u]llw]^' - «r^<+^) . 



2 



Now we introduce the entropy flux G\i'- 



r>n -pn+l ,,n+l J J + i ^+-1- J 



Then, by inserting (3.9) into (3.13) and by using the definition of G" i, one obtains 

' ' ' ' J+2 



u{v]n - u{v]) + Xi{g;^^ - g;_,) < -x,n] < o. 

This completes the proof of the proposition. D 



Next, we consider the entropy stability of the forward explicit Euler scheme (3.10). 



Hereafter, we focus on conservative schemes which admit the following viscosity form: 

("^) It"' ^ ^Hir^ - w. («'-^ft- - '='' - «'-^(- - ^-'') • 

so that the fiux /■ , i reads 

_ f{Vj+i) + f{Vj) 1 

h+l - 2 i^Qj+i{Zj+i -Zj). 

The matrix <5,+ i is a symmetric matrix whereas z = Uv{v) represent the entropy 
variables. It is easily seen that Z2 = Z2 = u and z^ = z^ = w. Here the conservative 
variables v are considered as functions of the entropy variables: in particular Vj = 
v{zj). The classical Lax-Friedrichs scheme and Rusanov scheme are particular cases 



of (3.15). Indeed, these schemes have the particular form 



(3.16) j^v, + — = — (P,,^{v,,,-v,)-p^_.{v,-v,^,, 
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with p,- 1 1 > a numerical viscosity. Indeed, one has 



Vzizj + ^{zj+i - Zj))dn {zj+i - Zj) = Vj+i 

Vz = {Zv)~'^ = iU^v)~'^ 



then Vz = vi. , so that (3.16) is a particular case of (3.15) by setting 



Qj+i = Pj^ I Vz{Zj + ^{Zj+i - Zj))d^. 



with Qj^i = Q'^,1- Following [Tj, one can compare the conservative scheme (3.15) 
with the entropy conservative scheme through the relation: 



"J'^i+t Jj-^/ -' i+l -^j-l 



(3.17) + - (^((%+i - %), D.^{Sj+i - %)) + {{zj - Sj_i), D^_i_{zj - Zj-_i))) 



with D-,1 = Q, I 1 — Q*,i where Q* i is the numerical viscosity matrix of the 



conservative scheme: 



1/2 
1/2 



2^C ( ^Z±i±^ + ^(5^.^, - 5,) ) d^ 



with C{z) = gz{z) and g{z) = f{v{z)) whereas J^j+i is a consistent numerical entropy 
flux. We prove the following proposition. 

Proposition 3.3. The (FE) finite difference scheme: 



:Q-+i(^7+i-^n- 



rn 



2-^i+iv J- 



zs entropy stable, i.e. there exists a numerical entropy flux Q"' i so that 

-'+2 

^«') - uiv]) + x^io;^, - g;_,) < o, 

under the following CFL condition 

m; (Aiiv;^i + A2||i?;+i||)' < Aimin(Sp(D;^p), 
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with M", A^", 1 defined as 



J + 5 



Mr 



^;^i 



sup ||f/™« + .« 



«G(0,1) 
-1 



n+1 
J 



-3))\\i 



\C{z-,^i{l-^,-~z-))\\di^\\Q 



i+ 



Proof. We first apply the Taylor Lagrange formula to U: 



n\T{„,n+l 



U{vJ+^) = U{v]) + U,{vjy {v] 



Jo 



Then, by using (3.17), one finds 



(3.19) - ^ (^1 - ^)^^;Vi(5^i - 5^) + (5; - 5}^i)^^^|(^ - 5;^i)) . 



Ai 

4 



The first term in the right hand side of (3.19) is positive and corresponds to entropy 



production due to the forward explicit Euler time discretization whereas the second 
term corresponds to entropy dissipation due to the spatial discretization. Next, we 
estimate the entropy production: in order to simplify notations, we set 



^ 



One has 



X"<- sup \\U,M + .{v' 



?6(0,1) 

Next, we estimate 11 f 



n+l 

j 



,n+l 



.ni|2 



m;ii<^,-. 



2' J ""3 



n||2 
j II • 



n+l 



f"|| by using (3.18): one finds 



|U,"+1 ,,"|| ^ \ II f" f" II I \ Ml R" II II ■y" ^"11 I 

II ^i - ^j II ^ ^lll/j+| - /j-l II + ^2 [W^j+l II Fi+l - ^i II + 

On the other hand, one has 



-D . 1 ^,' — Z 



•i-ii 



rn rn 



I' c [^ + m^i - ^)) d^ - q;+i) (5]Vi 



5;) 



^-5^i). 
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Then, by setting m^, = J^ \\C (5^" + ^{zj^^ - 5p) \\d^ + \\Q\i\\, one obtains 



II f" f" II <r Ar" II;;" ;:"ll m at'"- 11?" 5" II 
As a result, one finds that 

2 

II ~n _ ~n\\2 
ll^i+1 ^i II 



1^7 ^i-il 



27<M«(A,iv;^.+A,||5;^.| 

(3.20) +M" (AiAT" 1 + A2II5" 1 

Next, we set T" , = min (Sp(D"_^J). Furthermore, we assume that 

(3.21) M;(Aiiv;;_, + A2||5;+i||)' < XiV^^^.. 



Then, by using (3.21) together with (3.20) and (3.19), one obtains entropy stabihty 



for the exphcit forward Euler time discretization 
This completes the proof of the proposition. 



n 



Let us consider Lax-Friedrichs schemes: by applying proposition 3^ we prove the 
Corollary 3.4. Assume 

• there exists K > so that K~^ < M" < K for all j, n, 

• P;+i =p;_^i +max(|Sp(/.(t;;+i))|,|Sp(/,(t;p)|). 

6x 
The Lax Friedrichs type scheme, p" , = — -, is entropy stable under the CFL condi- 
tion 

for some constants Mj{K),j = 1,2. The Rusanov type scheme, p"" 1 = p > 0, is 
entropy stable under the CFL condition 

K{\^M^{K) + \2M2{K)f < pAi, 

Remark 3.5. The previous result states that the classical Lax-Friedrichs scheme is 
entropy stable if 6t = 0(5a:;^) whereas the Rusanov scheme is entropy stable only if 
6t = 0{5x^) which is the Von Neumann stability criterion found in section^ 
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4 Numerical Simulations 

4.1 Entropy stability: original vs new formulation 

Before carrying out a numerical simulation of a Liu-Gollub experiment with the full 



shallow water system (1.2), we have considered the more simple situation of a fluid 



over an horizontal plane without friction at the bottom. The shallow water system 
reads 



(4.1) 



dth + dx{hu) = 0, dt{hu) + dxQiu^ + g- 



f\j I V \y 'Y' ''p''p ' V % 



where g = 9.8m. s ^ and k = cr/p. The fluid under consideration in |LGj is an aqueous 
solution of glycerin with density p = 1.134 g.cm~^ and capillarity a = 67 dyn.cm~^. 



The augmented form of (4.1) reads 



(4.2) 



dth + dx{hu) = 0, 
dt{hu) + dx{hu^ + g 
dt{hw) + dx{huw) = 



h' 






We first tested the entropy stability of (second order accurate) difference approxima- 



tions for the shallow water equations (4.1) and for its "Schrodinger type" counterpart 
(4.2). We work on a finite interval of length X = 80cm with periodic boundary 
conditions. At time t = 0, the fluid velocity u = and the fluid height is given by 
h\t=o = h^ {l + 0.3exp(— 2000(x — 0.4)^)) with hjy = 1mm (the characteristic fluid 
height in Liu-Gollub experiments). In order to capture correctly the capillary ripples, 
we have chosen Sx = 0.25m,m. and St = 120Sx'^. In figure [T| we draw the profile of the 
surface of the fluid at time T = Is: as expected, the Lax-Friedrichs scheme slightly 
damps the capillary ripples in front of the shocks. 




Figure 1: Profile of the surface of the fluid at time T = 1 
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U\t 



as a function of time: the 



In figure 2, we have drawn the relative entropy 

picture clearly indicates that the difference approximation of (4.2 ) have better entropy 



stability properties than difference approximation of (4.1). 





Entro-pvas a functiDn uf time 












- origiinal fDrmulation: LF 

- original fo-rmulation: RU 

- new formulation: LF 

- new formulation: RU 


- 
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Figure 2: Entropy as a function of time: Comparison of the various spatial discretiza- 



tions for the new formulation of (4.1). As expected, the entropy is almost constant 



for the difference approximations of (4.2). The results for the Harten-Lax-van Leer 



scheme are not plotted as the results are the same than the ones for the Rusanov 
scheme 



A natural question arises about the new formulation: indeed one may ask whether 
the relation hw = ^y/^dxih^^"^) is satisfied for all time. If not, it does not make sense 
to compare the performance with respect to entropy stability since it would represent 
two distinct quantities. In figure [3} we draw the relative error at time T = 1 and 
defined as 



err^ 



\{hw) 



3/2 _ 3/2 
"j + 1 "j-1 

3Sx 



\hw\ 



J 



.N. 



The numerical simulations show very good agreement, especially for the less dissi- 
pative schemes, Rusanov and Harten-Lax-van Leer, than for Lax-Friedrichs scheme. 
Moreover, it is easily seen that at the linearized level, the difference approximations of 



(4.2) and of (4.1) have comparable CFL conditions: in particular the CFL condition 
for Rusanov scheme is of the form 5t = 0{5x'^^^) that is rather close to the "opti- 
mal" heuristic CFL condition 6t = 0((5x^). Therefore, we can conclude from these 



numerical simulations that a difference approximation of (4.2) with a Rusanov flux 
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and second order accurate both in time and space is a natural candidate to perform 
numerical simulations of Liu-Gollub experiments |LGj . 




Figure 3: Consistency of the new formulation: relative error between the new variable 
hw and ^y/ndx^h^^'^)- The Rusanov scheme and the Harten-Lax-Van Leer scheme 
have comparable consistency properties. In regards, the Lax Friedrichs scheme is less 
efficient to preserve consistency 



4.2 A semi-discrete entropy conservative scheme 

In this section, we use the Hamiltonian structure of the Euler-Korteweg equations to 
construct an entropy conservative scheme. For that purpose, we will write a semi dis- 
cretized form of the Euler Korteweg system which respects its Hamiltonian structure 
so that the entropy conservation is automatically satisfied. 



4.2.1 Derivation of a Hamiltonian difference approximation 

We consider the Euler-Korteweg equations with periodic boundary conditions and, 
for {g,u) = (pi, tij)j=i,...,Ar, we introduce the discrete Hamiltonian 

(4.3) H(,,u) = 5^p.| + F(p,) + ^/.(p„)(^^±i-^j . 



We also introduce the symmetric matrix J 

J = 



-In 
In 
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and the difference operator D, defined in the space of A^-periodic sequences in M^ as 
Dui = "'"'"2^"'"^ (the associated matrix D G M7v(M) is skew symmetric). Then, we 
introduce the Hamiltonian system 

d [ q\ ^ [ DVgH(^,u) 



More precisely, the discrete Hamiltonian system reads 

dpj {pu)j+i - ipu)j-i ^ 
dt 2Sx 

dUj ^ ^ ■+! - ^ ■-! ^ F'(p,+i) - i"'(p,-i) 



dt ASx 26 X 

+ 



1 ( f^'iPj+i) fpj+2 - Pj+iV i^'{pj-i) fpj - Pj-i^ ^ 




2 V 5x 

,Pj+i-pi 



-''^P^' 6x 

2..n— sx -^^^-^^^^^'^l-^- 



By construction, one has 



|h(,, u) = V,H(,, u)^^ + V,H(,, u)^^ 



H(f?, u) = V,H(f?, u)^ -^ + VuH(^, u) 

V,H(^, u)^ D V,H(f?, u) + V,H(f?, u)'^ D V,H(f?, u) 
V,H(^, u)^ (D + D^) VuH(^, u) = 0. 



It is clearly a consistent and first order discretization of the original system (1.3). At 
this stage, one can go one step further and derived naturally higher order entropy 
conservative scheme like in |CLj . Indeed, one easily improves the order of accuracy of 



(4.4) by considering a higher order approximation of the Hamiltonian and a higher 



order difference operator. 

As a consequence, one is left with the problem of finding a time discretization that 
preserves the Hamiltonian structure. It is easily seen that an explicit forward Euler 
time integration is unstable whereas the backward implicit Euler time integration is 
entropy stable. In the linearized case, the Crank Nicolson scheme preserves exactly 
the Hamiltonian. The hamiltonian difference scheme derived here is based on centered 
difference and it is well known that for hyperbolic conservation laws, this could be a 
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source of numerical instabilities or spurious oscillatory modes. Though, the scheme 
considered here also provide a control on the gradient of the density and thus on 
oscillatory modes in addition of being more stable. 



4.2.2 Numerical tests and dispersive shock waves 



In what follows, we have tested the difference hamiltonian approximation of (4.1). 



In order to be entropy stable it is necessary to employ an implicit method: we have 
used here an implicit backward Euler time discretization. Due to the nonlinearity 
of the problem, the Crank-Nicolson, second order accurate, time discretization does 
not guarantee entropy stability. Therefore, we did not try to compare with other 
schemes tested in the previous section. An important remark is that now there is 
no numerical viscosity: a drawback is that the dynamical behavior is dramatically 
changed as shown in figure 111 and [SJ) 



Formulalion namjitonjenne: hauleur defluide 





Figure 4: On top: Fluid height profile at time T = 0.26: after the breaking of 
the wave, an oscillatory part appears, due to dispersive regularization. On below: 
Fluid height profile at time T=0.66: the oscillatory zone is wider, characteristic of a 
dispersive shock wave 

In order to see whether it is a numerical artifact, we checked the entropy stability of 
the difference hamiltonian approximation: the entropy remains clearly bounded with 
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time. Indeed, these oscillations are not a numerical artifact and can be explained 
(formally) by the theory of dispersive shock waves. Here, the classical hyperbolic 
shocks are smoothed by disperses effects: the oscillatory zone grows up in time and 
the oscillations are described by the Whitham modulations equations. This picture 
is not valid anymore in the presence of a slight amount of viscosity: there are still 
some oscillations but the width of the oscillatory zone stops growing after some time: 
see |J] and |EGKj for a detailed analysis respectively in the case of the Korteweg de 
Vries-Burgers equation and in the case of the Kaup system perturbed by a viscous 
term. In the previous section, a similar situation arises: the oscillatory zone stops 
growing as it is shown in the former computations in picture [l} Here the physical 
viscosity is replaced by numerical viscosity. 




Figure 5: Entropy as a function of time. The behavior is rather similar to difference 
approximations in original variables: the entropy first increases then decreases with 
time: the "holes" in the decreasing part of the curve correspond to times when the 
bumps interact 



4.3 Simulation of a Liu-Gollub experiment 

In this section, we show a numerical simulation for a shallow water model derived 
for thin film flows down an inclined plane. The model lies in the family of first order 
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models derived in |NV] : it is written as 
(4.5) 

dth + dx{hu) = 0, 

dt{hu) + d^ [hu^ + P(/i, Ai)) = Ai (ghsm{d) " 3^ + -hd^^^h) + Aud^^q, 

where Ai > and P{h, Ai) is a pressure term given by 

r^/T . \ . //ix^^ /4 2y4i, /5fsin6'\ , . 

P(;,^)=AscoB(9)- + (--— )(^^jft=. 

Note that the viscous term AudxxQ is only heuristic so that the model is not a second 
order accurate model (with respect to the aspect ratio e = H/X -C 1, if pa 1mm, A ~ 
1cm). Here, the viscosity plays a significant role, z/ = 6.28 cS*. Indeed, the Reynolds 
number for this experiment is rather low: Re = 29 whereas the slope of the bottom 
is set to 6 = 6.4°. The frequency of the perturbation at the inlet is / = 1.5Hz. 
We have chosen to carry out the numerical simulations with Ai = 1. Following the 
conclusions of our study, we have chosed to carry out numerical simulations with a 



fully second order accurate scheme of the extended formulation of (4.5) and used a 
Rusanov fiux for the first order part. Our numerical results show a good agreement 
with the experiment by Liu and GoUub |LGj . 

Up to now, the choice of boundary conditions for the Euler-Korteweg equations on a 
finite interval is an open problem so that we have chosed rather arbitrary boundary 
conditions. Furthermore, since the difference scheme contains numerical/physical 
viscosity, we have considered a set of 5 boundary conditions. First, at the inlet, we 
chosed: 

h\x=o = hN {I + 0.03 sm{27!- ft)) , q\x=o = qN, dxh\x=o = 0. 

In contrast to |KRS Vj . we have chosed free boundary conditions at the outlet: 

dxh\x=L = dxq\x=L = 

instead of "hyperbolic type" boundary conditions where h and q are advected with 
an artificial velocity Vout > 0. As pointed out in |KRSVj the choice of the boundary 
conditions at the outlet does not seem to influence the dynamic within the channel 
(no reflection waves). 



5 Concluding remarks 

In this paper, we considered the stability of various difference approximations of the 
Euler-Korteweg equations with apphcations to shallow water equations with surface 
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Figure 6: Simulation of Liu Gollub experiment |LGj : the Reynolds number is Re = 
29 and the inclination is 6* = 6.4°. The frequency at the inlet if / = 1.5Hz. On top: 
a picture of the complete experiment, from the inlet to the outlet (2m). On below: a 
zoom over one spatial period when roll-wave profiles are stabilized 



tension. A first class of difference approximations is built by considering the Euler- 
Korteweg system as the classic isentropic compressible Euler equations perturbed by 
a disperse term. This latter term is discretized with centered finite differences and 
various classical scheme for the convection part are considered. It is proved that a 
certain amount of numerical viscosity is needed to obtain difference schemes that are 
stable in the von Neumann sense (under suitable CFL conditions). 

In order to get entropy stability, we considered an extended formulation of the Eu- 
ler Korteweg equations and proved entropy stability of Lax-Friedrichs type schemes. 
We have shown numerically that the extended formulation of the Euler-Korteweg 
system has better stability properties than the original one. We also carry out a 
numerical simulation of a shallow water system which models an experiment by Liu 
and Gollub to observe roll- waves |LGj . 

By considering the Euler-Korteweg system as a Hamiltonian system of evolution 
PDEs, we introduced a semi-discretized difference approximations which preserves 
the Hamiltonian structure. This scheme has no numerical viscosity so that it is 
particularly useful to study purely dispersive Euler Korteweg system: in particular. 
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one can find numerically the dispersive shock waves [E] of the Euler-Korteweg system. 

Several questions remain open. First, we carried out a numerical simulation of 
the Liu GoUub experiment by choosing rather arbitrary boundary conditions. In 
fact, the choice of suitable boundary conditions for the Euler Korteweg system on 
a finite interval in order to prove well posedness is still an open problem. A first 
attempt in this direction is found in |X] where the well posedness of the linearized 
Euler-Korteweg equations is proved on a half space under a generalized Lopatinskii 
condition. 

Furthermore, we restricted our attention to one dimensional problem. For thin 
film flows, this restricts the study to primary instabilities: in order to analyze sec- 
ondary instabilities found experimentally |LSG] . one has consider 2d problems. In 
that setting, an extended formulation is still available |BDDdj so that we expect our 
analysis extends easily, at least to cartesian meshes. An other interesting question is 
the extension of this analysis to other mixed hyperbolic/dispersive equations like the 
Boussinesq equations or the Serre/Green-Naghdi equations. Up to now, the strategy 
adopted to deal with these system is time splitting without proof of stability (though 
numerical results are rather satisfying). 

Finally an other open interesting question concerns the time integration of the 
hamiltonian semi-discrete approximation: here, we have used a backward implicit 
Euler time integration so as to be entropy stable but it does not preserve the hamil- 
tonanian (nor a perturbation of it). This kind of method are particularly of interest 
in order to study the nonlinear stability of various traveling waves solutions of the 
purely dispersive equations. 

References 

[A] C. Audiard Non homogeneous boundary value problems for linear dispersive 

equations, Comm. Parti. Diff. Eqs 37 (2012) no 1, p. 1-37 

[BDD] S. Benzoni-Gavage, R. Danchin, S. Descombes Well-posedness of one- 
dimensional Korteweg models, Electronic J. Diff. Equations 59 (2006), pp. 
1-35. 

[BDDd] S. Benzoni-Gavage, R. Danchin, S. Descombes On the well-posedness for 
the Euler-Korteweg model in several space dimensions, Indiana Univ Math 
Journal 56 (2007), no 4, pp. 1499-1579. 



REFERENCES 33 

[CL] C. Chalons, P.G. LeFloch High- Order Entropy- Conservative Scheme and Ki- 
netic Relations for Van der Waals Fluids, Journal of Computational Physics 
168 (2001) 184-206. 

[Col] F. Coquel, P.G. LeFloch An entropy satisfying MUSCL scheme for systems 
of conservation laws Numer. Math. 74 (1996) 1-33. 

[E] G. A. El, Resolution of a shock in hyperbolic systems modified by weak dis- 

persion Chaos, 15 (3): 037103, 21, 2005. 

[EGK] G. A. El, R. H. J. Grimshaw, and A.M. Kamchatnov, Analytic model for a 
weakly dissipative shallow water undular bore Chaos, 15 (3): 037102, 2005. 

[EGS] G. A. El, R. H. J. Grimshaw, and N. F. Smyth. Unsteady undular bores 
in fully nonlinear shallow-water theory. Physics of Fluids, 18(2):027104, 
February 2006. 

[J] S. Johnson A non-linear equation incorporating damping and dispersion J. 

Fluid. Mech. 42 (1970) p. 49-60. 

[KRSV] S. Kalliadasis, C. Ruyer-Quil, B. Scheid, M.G. Velarde Falling Liquid Films, 
Applied Mathematical Sciences 176 (2012) 

[LMR] P.G. LeFloch, J.M. Mercier, C. Rohde Fully discrete, entropy conservative 
schemes of arbitrary order, SIAM J. Numer. Anal. 40 (2002) No 5, pp 1968- 
1992. 

[LG] J. Liu, J. P. Gollub Solitary wave dynamics of film flows, Phys. Fluids 6, 1702 
(1994). 

[LSG] J. Liu, B. Schneider, J. P. Gollub Three-dimensional instabilities of film flows, 
Phys. Fluids 7, 55 (1995). 

[NV] P. Noble, J. -P. Vila Consistency, Stability, Callilean indifference criteria for 
design of Two moments closure equations of Shallow water type for Thin 
Film laminar flow gravity driven, in preparation. 

[VL] B. Van Leer Towards the ultimate conservative difference schemes: a second 
order sequel to Codunov's method J. Comp. Phys 43 (1981) 357-372. 

[T] E. Tadmor Entropy stability theory for difference approximations of nonlin- 

ear conservation laws and related time- dependent problems Acta Numerica 
(2003) pp. 451-512. 



